A fusion of the LAPW and the LMTO methods: the augmented plane wave plus 

muffin-tin orbital (PMT) method 
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We present a new full-potential method to solve the one-body problem, for example, in the local 
density approximation. The method uses the augmented plane waves (APWs) and the generalized 
muffin-tin orbitals (MTOs) together as basis sets to represent the eigenfunctions. Since the MTOs 
can efficiently describe localized orbitals, e.g, transition metal 3d orbitals, the total energy con- 
vergence with basis size is drastically improved in comparison with the linearized APW method. 
Required parameters to specify MTOs are given by atomic calculations in advance. Thus the robust- 
ness, reliability, easy-of-use, and efficiency at this method can be superior to the linearized APW 
and MTO methods. We show how it works in typical examples, Cu, Fe, Li, SrTiOs, and GaAs. 

PACS numbers: 71.15.Ap, 71.15.Fv 71.15.-m 



There are several kinds of all-electron full poten- 
tial (FP) methods to obtain numerically-accurate so- 
lutions in the local density approximation to the den- 
sity functional theory Among such FP methods, 
most popular ones are the linearized augmented plane 
wave (LAPW) method, and the projector augmented- 
wave (PAW) method 0, S H H- They both use plane 
waves (PWs) to expand the eigenfunctions in the inter- 
stitial regions. However, PWs do not efficiently describe 
the localized character of eigenfunctions just around the 
muffin-tins (MTs). For example, oxygen 2p (denoted as 
0(2p) below) and transition metal's 3d orbitals are well 
localized and atomic-like even in solids, thus we need to 
superpose many PWs to represent these orbitals. For ex- 
ample, as shown in Refla (and below), the energy cutoff 
for the augmented PW(APW) > 15Ry is required 
in fee Cu to obtain ~ 1 mRy convergence for total energy 
in LAPW. In contrast, such orbitals can be rather effec- 
tively represented by localized basis in real space. In fact, 
it is already accomplished in the linearized muffin-tin or- 
bital (LMTO) method, which differs from the LAPW 
method in that envelope functions consists of atomic-like 
localized orbitals @, Hf instead of PWs. Such a localized 
augmented waves are called as MT orbital (MTO). 

To circumvent the inefficiency in the LAPW method, 
we have implemented a new method named as linearized 
augmented plane wave plus muffin-tin orbital (PMT) 
method. The PMT method includes not only the APWs 
but also MTOs in its basis set. Our implementation be- 
comes LAPW in the no MTO limit. As we show later, 
we can very effectively reduce the number of basis set 
by including MTOs; we see the rapid convergence of the 
total energy as a function of the number of APWs (or 
energy cutoff -Emax)- As f ar as we tested, APWs with 
-EjmAX ~ 5 Ry in addition to minimum MTOs will be rea- 
sonably good enough for usual applications; e.g. for <1 
mRy convergence of total energy for Cu. Even in compar- 
ison with the LMTO method of Ref.@], the PMT method 
is quite advantageous in its simplicity. The parameters 
to specify these minimum MTOs (-Eh and Rn for each I 
channel. See next paragraph.) are automatically deter- 



mined by atomic calculations in advance. This is a great 
advantage practically because optimization of these pa- 
rameters is a highly non-linear problem Q which makes 
the LMTO difficult to use. Thus the PMT method can 
satisfy the requirements for latest first-principle meth- 
ods, reliability, speed, easy-of-use, and robustness very 
well. In what follows, we explain points in our method, 
and then we show how it works. 

We adapted a variant of the LMTO method developed 
by Methfessel, van Schilfgaarde, and their collaborators 
0, Q . This method uses the "smooth Hankel function" 
invented by Methfessel as a modification of the usual 
Hankel functions so as to mimic the atomic orbitals [8|, l9( . 
It contains an additional parameter called as the smooth- 
ing radius i?H in addition to the usual energy parameter 
En which specifies the damping behavior. By choosing 
i?H, we can control bending of the function just outside 
of MT. By using this degree of freedom, we can reproduce 
eigenvalues of atoms very well even if we substitute the 
eigenfunction outside of MT with such a smooth Han- 
kel function. This is an important feature of the func- 
tion in comparison with others like Gaussians, which are 
not directly fit to the atomic orbitals. Analytic proper- 
ties of the smooth Hankel are analyzed in detail by Bott 
et al, and all the required operations, e.g, Bloch sum, 
to perform electronic structure calculations are well es- 
tablished [§]. The augmentation procedure requires the 
one-center expansion of the envelops functions. In our 
method, the one-center expansion is given as the linear 
combinations of the generalized Gaussian (Gaussians x 
polynomials), where the expansion coefficients are deter- 
mined by a projection procedure as described in Sec. XII 
in Ref. 9. Then the generalized Gaussians in each angular 
momentum /m-channel are replaced by the linear combi- 
nations of a radial function (pi and its energy derivative 
4>i in the usual manner of augmentation |8|. When we 
use high energy cutoff -E'max' ~ 15Ry, we needed to use 
~ 15 generalized Gaussians for the one-center expansion; 
however, ~ 5 is good enough for practical usage with 

Another key point in our method is the smoothed 
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FIG. 1: (color online) Total energy is plotted as a function of (number of basis at F point -/Vj) -1 - Lower panels for each 
material zoom up upper panels in energy (in SrTiOa, data points in LAPW are out of the scale in lower panel). We also plot 
spin moment (/xb) for Fe, and the band gaps for GaAs and SrTi03 by dotted lines, referring to right y-axis. The absolute 
values of energies are ambiguous (See text in Table [TJ) . -Emax(R-y) are shown next to each data point. "Ad lo" means we treat 
4d as local orbital. PMT denotes our new method with APW+MTO, where we use minimum MTO basis; number of MTOs 
(including lo) are 14(Cu),5(Li),9(Fe),18(GaAs), and 33(SrTi0 3 ). Number of k point in 1st Brillouin zone (BZ) are 12 3 for Cu 
and Li, 20 3 for Fe, and 6 3 for GaAs and SrTiOa. Lattice constants are not necessarily at their total energy minima (e.g., 6.8a.u 
for Cu). See Table U for exchange-correlation used. 



charge density treatment introduced by Soler and 
Williams [H to the LAPW method. The charge den- 
sity is divided into the smooth part, the true density on 
MTs, and the counterpart on MTs. The smooth part 
is tabulated on regular uniform mesh in real space, and 
the others are tabulated on radial mesh in the spheri- 
cal harmonics expansion in each MT sphere. The PAW 
methods [1, [f| also use this treatment. It enables low- 
angular momentum I cutoff for augmentation and makes 
the calculation very efficient. As for the regular mesh in 
real space, we usually need to use the spatial resolution 
corresponding to the cutoff energy -EJ^ax* 1 = 10 ~ 15 Ry, 
which is determined so as to reproduce the smooth Han- 
kel function well in real space. Note that we still have 
some inefficiency, e.g. an atom in a large supercell re- 
quires a fine mesh everywhere only in order to describe 
the density around the atom. This problem is common 
to any method which uses an uniform mesh for density. 

Though the LMTO formalism shown in 0, H is in- 
tended for such MTOs constructed from the smooth Han- 
kel functions, it is essentially applicable to any kind of en- 
velope functions. As we explained above, our formalism 
is not so different from LAPW/PAW formalisms shown 
in 0, [1[ except in the augmentation(projection) proce- 



dure. The atomic forces are calculated [8|,ll0( in the same 
manner as in PAW [5j. For deep cores, we usually use a 
frozen core approximation which allows the extension of 
core densities outside of MT (but no augmentation) 0]. 
Further, we use the local orbital (lo) method [ll| in some 
cases; for example, to treat 3d semicore for Ga (denoted 
as Ga(3<i[lo]) below), or to reproduce high- lying bands 
for Cu by Cu(4d[lo]). 

Results: In Fig[TJ we plot the total energies as func- 
tions of the inverse of the number of basis at the T point 
(7V^) _1 in order to observe its convergence as N£ — ► oo 
(the number of basis is controlled by -Emax as shown on 
Fig[T]together). Here we includes minimum MTOs whose 
parameters En and i?H are fixed by the atomic calcula- 
tions (here "minimum" means only the atomic bound 
states). The lo's are included as explained in Fig[T] 

There is a problem of linear dependency in the basis 
set when we use large E^ax- For example, in Li, we 
could not include MTOs of Li(2s2p) (Li(2s) and Li(2p)) 
as basis at E^^[ > 6, because then they are well ex- 
panded by PWs. This occurs also for other cases; when 
we use B^AX high enough to expand a MTO, the rank 
of the overlap matrix of basis set is reduced by one. Thus 
possible ways to use large i?^AX are: (1) keep only well 
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localized MTOs which are not yet expanded by given 
^MAJC' or (2) remove a subspace of basis through the 
diagonalization procedure of the overlap matrix before 
solving the secular equation. For (2), we need to intro- 
duce some threshold to judge the linear dependency. This 
can cause an artificial discontinuity when changing lat- 
tice constants and so on. Thus (1) should be safer, but 
here we use the procedure (2) with careful check so that 
such discontinuity do not occur. We use the number of 
basis after the procedure (2) for N£ to plot Fig. [TJ Even 
for LAPW cases, we applied the procedure (2), e.g, to 
the case for £"max =20Ry in SrTiOs; then we reduce the 
dimension of Hamiltonian from 606 to Ng = 550 (data 
point at the left end of SrTi0 3 in FigHJ. 

As is seen in all the cases, the PMT method shows 
smooth and rapid convergence for the total energy at 
(iV^) _1 — > 0. On the other hand, the convergence in 
LAPW (no MTO limit in our PMT implementation) 
shows a characteristic feature; it is way off until it reaches 
to some i?MAX> e -S> ~ ^ tyf m Cu. This ^ s consistent 
with previous calculations p. ^majc ^ 15 is required 
to reproduce the 3d localized orbitals. We also show the 
convergence behaviors for band gap and magnetic mo- 
ments by dotted lines (right y-axis); they are quite satis- 
factory. 




FIG. 2: Energy bands for Cu up to 100 eV above the Fermi 
energy. Top-left is for APW only without local orbitals. 
Top-right: APW +4d[lo]. Bottom-left: APW plus MTO (9 
basis) +4d[lo]. Bottom-right: MTO(34 basis) +4d[lo]. All are 
converged for -Emax- 

Figf2] shows the energy bands up to lOOeV above the 
Fermi energy. Though the Cu panel in Figfl] shows the 
little effects of Ad local orbital to the total energy, it 
affects energy bands above ~ 30eV. Calculations includ- 
ing 4d[lo] gives good agreements each other. This means 
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FIG. 3: Total energy v.s. lattice constant. Labels "Emax" 
means 

#max ■ MTO (high) for Cu is explained in text. 



that we have no artificial bands (ghost bands). The 
MTO (high) panel is by the pure LMTO method where 
we use 34+5 basis (spdfg + spd+Ad[lo\). With some care- 
ful choice of Er and i?H , the LMTO method can be very 
efficient and accurate. 

n(MTO) 
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n(APW) 

FIG. 4: Total energy for different type of MTO with chang- 
ing the number of APWs n(APW) for SrTi0 3 . Each 
curves(A,B,C,D,E,F) corresponds to a type of MTO sets (see 
texts). For example, curve C uses 36 MTOs (including lo- 
cal orbital). Thus curve C at n(APW)=80 means a calcu- 
lation with n(MTO)+n(APW)=36+40=76 basis. The con- 
verged energy ~ —2.76 Ry is a little different from the Fig[T] 
because this calculation is performed with coarse conditions 
for test purpose. 

Fig[3] shows the total energy as function of the lattice 
constant. All lines looks very parallel. This shows that 
PMT do not include some systematic errors. This is true 
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TABLE I: Calculated lattice constant, bulk modulus, and 
cohesive energies. Values are by the PMT method with 
#ma¥ = 15Ry (£ma¥ = 6Ry for Li). Values in parenthe- 
sis are with -Bmajc = 5Ry. We show values from other lit- 
eratures together. We used the Barth-Hedin exchange cor- 
relation [3] except Fe where we use the VWN type 
Since we have not determined the total energies of each atoms 
exactly (we assume spherical atoms and assume some elec- 
tron configurations), the cohesive energies (constant in en- 
ergy axis in Fig. Q] and Fig. [3]) are somehow ambiguous; espe- 
cially for cases including transition metals. For comparison 
with other calculations, it will be better to use total energies; 
add -3304.4345Ry(Cu), -14.7106Ry(Li) ,-2540.4767Ry(Fe), - 
8397.5970Ry(GaAs), and -8502.4637Ry(SrTiO 3 ) to the cohe- 
sive energies. 





Lattice constant. 


Bulk Modulus 


Cohesive 




(a.u.) 


(GPa) 


energy(Ry) 


fee Cu 


6.650(6.649) 


188(187) 


-.332(-.331) 


LAPW 


6.65 


192 




expt. a 


6.81 


131 




bec Fe 


5.209(5.208) 


258(259) 


-.667(-.665) 


PAW 6 


5.20 


247 




LAPW 


5.210 


245 




expt. c 


5.417 


172 




bec Li 


6.347(6.341) 


15.3(15.5) 


-.124(-.123) 9 


PAW 6 


6.355 


15.0 


-.149 


SrTi0 3 


7.367(7.360) 


220(225) 


-2.731 (-2.709) 


PP d 


7.31 


203 




expt. e 


7.39 


184 




GaAs 


10.61(10.61) 


74.9(75.0) 


-.576(-.574) 


LAPW 


10.62 


74 


-.587 


expt. e 


10.68 


76 


-.479 



"Ref. 
6 Ref. 
c Ref. 



c Ref. 
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''Ref. 16 | . PP denotes a pseudopotential method. 



9 The difference from PAW may be because it uses the non- 
polarized Li atom as reference. If we do so, we have -.150(-.150). 



also for other materials (not shown) . Table [J shows the 
obtained lattice constants and related parameters. We 
also showed values for E^f^— 5Ry. Even though we 
still need other extensive tests, we think that such low 
cutoff is reliable enough for practical applications. For 
^max = 5Ry, we only need ~ 25 basis (when we do not 
include lo) for Cu. 

FigfJ] shows the total energies with different MTO sets 
combined with different numbers of APWs. Curve "A" 
includes just MTOs for 0(2s2p), sufficient for a crude 
representation of the valence bands. Also included are 
Sr(4p[lo]) and Ti(3p[lo]). A large number of APWs is 
needed to get a good total energy: ~ 150 APWs are 
needed to converge energy to within ^50 mRy of the 
converged result. "B" corresponds to an extreme tight- 
binding basis, consisting of Sr(5s5p) and Ti(4s4p4d) in 
addition to "A" . (Note that the conduction band is 
mainly Ti(3d), and 0(2s2f>)). The total energy of the 



MTO basis alone (no APWs) is rather crude — more 
than 200 mRy underbound. However, only 25 orbitals 
(plus 6 for lo's) are included in this basis. The energy 
drops rapidly as low-energy APWs are included: adding 
^40 APWs is sufficient to converge energy to ~50 mRy. 
As more APWs are added, the gain in energy becomes 
more gradual; convergence is very slow for large -E'max'' 
"C" differs from "B" only in that Sr(3d) orbital was 
added. With the addition of these 5 orbitals, the MTO- 
only basis is already rather reasonable. This would be 
the smallest acceptable MTO-basis. As in the "B" basis, 
there is initially a rapid gain in energy as the first few 
APWs are added, followed by a progressively slower gain 
in energy as more APWs are added. "D" is a standard 
LMTO minimum basis: spd orbitals on all atoms. Com- 
paring "C" or "E" to "A" shows that the MTO basis is 
vastly more efficient than the APW basis in converging 
the total energy. This is true until a minimum basis is 
reached. Beyond this point, the gain APWs and more 
MTOs improve the total energy with approximately the 
same efficiency, as the next tests show. "E" is a stan- 
dard LMTO larger basis: spd + spd orbitals on Sr and 
Ti, and spd + sp on O. Comparing "D" and 'F" shows 
that the efficiency of any one orbital added to to the 
standard MTO minimum basis is rather similar in the 
APW and MTO cases. Thus, increasing the MTO basis 
from 51 to 81 orbitals in the MTO basis lowers the en- 
ergy by 33 mRy; adding 33 APWs to the minimum basis 
("D") lowers the energy by 36 mRy. "F" enlarges the 
MTO basis still more, with Sr: spd + spd, Ti: spd + spd, 
O: spd + spd. Also local orbitals are used to represent 
the high- lying states Ti(4d[lo]) and 0(3s[lo],3p[lo]). For 
occupied states, these orbitals have little effect for total 
energy as in the case of Cu. 

In conclusion, we have implemented the PMT method 
whose basis set consists of the APWs together with the 
MTOs which are localized in real space. This idea is con- 
sistent with the nature of the eigenfunctions in solids; 
they can be localized as atoms or extended as PW. 
This method combines advantages of LAPW/PAW and 
LMTO. We have implemented force calculations, but no 
results are shown here. One of the advantage is in its 
flexibility. As an example, in order to treat Cu impu- 
rity in Si, it will be possible to choose very low SjyfAx 
just responsible for empty regions because MTOs for Cu 
and MTOs for Si span its neighbors already very well. 
Convergence is easily checked by changing £"max (much 
simpler than LMTO). In future, we can use the PMT 
method to make a natural division of the Kohn-Sham 
Hamiltonian into localized blocks and extended blocks, 
instead of the energy windows method for the maximally 
localized Wanner functions [lg|. The problem of large 
E UAX should be solved to make PMT more efficient. 

This work was supported by ONR contract N00014-7- 
1-0479. We are also indebted to the Ira A. Fulton High 
Performance Computing Initiative. 



5 



[1] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 
[2] O. K. Andersen, Phys. Rev. B 12, 3060 (1975). 
[3] J. M. Soler and A. R. Williams, Phys. Rev. B 40, 1560 
(1989). 

[4] P. E. Blochl, Phys. Rev. B 50, 17953 (1994). 

[5] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 

[6] A. Y. Liu, D. J. Singh, and H. Krakauer, Phys. Rev. B 

49, 17424 (1994). 
[7] M. Mesfessel and M. van Schilfgaarde, 'NFP manual 1.01 

Oct 10,1997'. NFP is previous to the current LMTO 

package lmf maintained by M. van Schilfgaarde. 
[8] M. Methfessel, M. van Schilfgaarde, and R. A. Casali, in 

Lecture Notes in Physics, edited by H. Dreysse (Springer- 

Verlag, Berlin, 2000), vol. 535. 
[9] E. Bott, M. Methfessel, W. Krabs, and P. C. Schmidt, J. 

Math. Phys. 39, 3393 (1998). 
[10] M. Methfessel and M. van Schilfgaarde, Phys. Rev. B 48, 

4937 (1993). 



[11] D. J. S. E. Sjostedt, L. Nordstrom, Solid State Commu- 
nications 114, 15 (2000). 

[12] U. von Barth and L. Hedin, J. Phys. C 5, 1692 (1972). 

[13] L. W. S. H. Vosko and M. Nusair, Can. J. Phys. 58, 1200 
(1980). 

[14] A. Khein, D. J. Singh, and C. J. Umrigar, Phys. Rev. B 
51, 4105 (1995). 

[15] L. Stixrude, R. E. Cohen, and D. J. Singh, Phys. Rev. B 
50, 6442 (1994). 

[16] L. M. Liborio, C. G. Sanchez, A. T. Pax- 
ton, and M. W. Finnis, Journal of Physics: 
Condensed Matter 17, L223 (200 5), URL 
|http : //stacks . iopTor g/0953- 8984/ 17/L223| 

[17] C. Filippi, D. J. Singh, and C. J. Umrigar, Phys. Rev. B 
50, 14947 (1994). 

[18] I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 
65, 035109 (2001). 



